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Abstract. We present analytical and numerical studies of magnetorotational instabilities occuring in magnetized 
accretion disks. These calculations are performed for general radially stratified disks in the cylindrical limit. We 
elaborate on earlier analytical results and confirm and expand them with numerical computations of unstable 
eigenmodes of the full set of linearised compressible MHD equations. We compare these solutions with those 
found from approximate local dispersion equations from WKB analysis. 

In particular, we investigate the influence of a nonvanishing toroidal magnetic field component on the growth rate 
and oscillation frequency of magnetorotational instabilities in Keplerian disks. These calculations are performed 
for a constant axial magnetic field strength. We find the persistence of these instabilities in accretion disks close 
to equipartition. Our calculations show that these eigenmodes become overstable (complex eigenvalue), due to 
the presence of a toroidal magnetic field component, while their growth rate reduces slightly. 
Furthermore, we demonstrate the presence of magneto-rotational overstabilities in weakly magnetized sub- 
Keplerian rotating disks. We show that the growth rate scales with the rotation frequency of the disk. These 
eigenmodes also have a nonzero oscillation frequency, due to the presence of the dominant toroidal magnetic field 
component. The overstable character of the MRI increases as the rotation frequency of the disk decreases. 
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1. Introduction 

Accretion disks are a rather common phenomenon in as- 
trophysics. They vary in size from a few up to a hundred 
astronomical units (AU) in young stellar objects (YSO) 
to approximately a hundred parsecs (pc) in the centers of 
active galactic nuclei (AGN). The central accreting object 
is a protostar in YSOs, a white dwarf, neutron star or a 
black hole in a binary, and a massive black hole in the cen- 
ter of an AGN. Part of the gravitational energy released 
during the accretion process is radiated away over a wide 
range of frequencies. Multi-wavelength studies obtained 
from modern observational facilities have significantly in- 
creased our knowledge about the nature of accretion disks 
and their associated central objects. 

On the other hand, insights into the nature of the ac- 
cretion process itself have largely been obtained by theo- 
retical and computational studies. In all the above men- 
tioned astrophysical objects, one needs outward transport 
of angular momentum of the accreting material in order 
to sustain an accretion disk around the central object. 
This angular momentum transport can be sustained by a 
turbulent viscosity mechanism operating inside the disk 
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material itself, where this mechanism excerts a torque on 
the accretion disk (Shakura & Sunvaev 11973}! . This tur- 
bulent viscosity, in turn, can originate from the devel- 
opment of fluid or magnetofluid instabilities occuring in 
the accretion disk. In the early 1990s it was realised by 
Balbus & Hawley ( 199*U that magneto-rotational instabil- 
ity (MRI) could provide the physical basis for this angular 
momentum transport in accretion disks. This instability 
was already known in the literature (Velikhov 1959 and 
Chandrasekhar 1960), but had not been applied in the 
context of accretion disks. 

In more recent years, global magnetohydrodynami- 
cal (MHD) simulations of accretion disks have been per- 
formed, where much of the dynamics is interpreted as a di- 
rect consequence of the presence of the magneto-rotational 
instability (see for example Hawley, Balbus & Stone 
2001). Other papers rather claim that both convective 
and magneto-rotational instabilities can play a dominant 
role in the transport of angular momentum (see for exam- 
ple Igumenshchev, Narayan & Abramowicz 2003 ). These 
types of disks are referred to as convection-dominated ac- 
cretion flows (CDAF). Less attention has been paid in 
recent years to the spectral analysis of instabilities oc- 
curring in accretion disks (see however Christodoulou, 
Contopoulos & Kazanas 2003). 



Ihe aim or this paper is to present a detailed linear 
analysis of magnetorotational instabilities present in a va- 
riety of global accretion disk configurations. We will con- 
sider a sample of accretion disk configurations that vary 
from sub-Keplerian (thick) to Keplerian (thin) rotating 
disks. In the case of Keplerian rotating disks, we will con- 
sider both weakly magnetized disks and disks that are 
close to equipartition. In particular, we investigate the 
influence of a nonvanishing toroidal magnetic field com- 
ponent on the growth rate and oscillation frequency of 
magnetorotational instabilities in Keplerian disks. To our 
knowledge, this has not yet been done. 

In this paper, we will limit ourselves to axisymmet- 
ric instabilities and the configurations are all taken in 
the cylindrical limit. The calculations continue the work 
presented by Keppens, Casse & Goedbloed (2002), ad- 
vocating the need for a more elaborate magnetohydro- 
dynamic spectroscopic analysis of all waves and instabil- 
ities in magnetized disks. We will use a semi-analytical 
approach, as well as numerical solutions to the full set of 
linearised, compressible MHD equations obtained with the 
code LEDAFLOW (Nijboer et al. EWI) . 

We will first show that magneto-rotational instabilities 
are present in both sub-Keplerian and Keplerian rotating 
disks and compare the growth rates obtained for these 
models. We find the presence of MRI even for cases where 
the disk is close to equipartition (i.e. /3 ~ 1), however 
the strength of the toroidal magnetic field should then be 
much larger than the axial magnetic field strength. 

This paper is organised as follows: in section 2, we 
recall the essential elements from spectral theory of MHD 
waves and instabilites. In section 3, we present the model 
and the limitations of the accretion disk configuration we 
use. In section 4, we explain the numerical strategies. In 
section 5, we discuss the magnetorotational overstabilitics 
which occur in our models and finally, in section 6, we 
summarise and present our conclusions. 

2. Spectral theory 

We make use of the ideal MHD approximation to model 
an accretion disk in the presence of a magnetic field. This 
assumes that the disk matter consists of a plasma that is 
sufficiently ionized and that one treats the dynamics on a 
length scale such that the one-fluid approximation is valid. 
The ideal MHD equations are 

p— = -pvVv- Vp+j xB-pV$, (1) 
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dB 
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where the variables, p, p, v, B, j = V x B and 7 
are the density, pressure, velocity, magnetic field, gravi- 
tional potential, electric current and ratio of the specific 



heats, respectively, furthermore, from Maxwell s theory, 
the equation V • B = must also be satisfied. We have 
assumed a dimensionalization where the permeability of 
vacuum po = 1. Equations Q-® are the momentum, en- 
tropy, induction and mass conservation equation, respec- 
tively. 

2.1. Frieman-Rotenberg formalism 

To investigate the stability properties of accretion disks, 
we linearize the ideal MHD equations ©-©. The lin- 
earisation is done by assuming a time-independent equi- 
librium. This assumption is justified at time scales which 
are much shorter than the dynamical time scale of the 
disk equilibrium. When linearizing, we write all variables 
of the ideal MHD equations in the following way: 



/ = /o(r) + /i(r,t), and|/ 1 |«/ . 



(5) 



Here, /o is the equilibrium quantity while f± is the time- 
dependent fluctuation about the equilibrium quantity. 
The resulting equations can be rewritten in terms of the 
Lagrangian displacement field, £(r, t), which is related to 
the perturbed velocity, vi, in the following way 



v ' = ( §1 + v o • v U - i ■ Vv <>- 



(6) 



where Vo represents the equilibrium flow velocity. By in- 
troducing the Lagrangian displacement field, there is no 
longer confusion between equilibrium and perturbed quan- 
tities when we suppress the subscript on all equilibrium 
quantities. The govering equation for the displacement 
field, called the Frieman-Rotenberg equation (Frieman 
and Rotenberg 11960(1 . is 



P0 + ^.V§-F«)= Ql 

where F(£) represents the force operator, defined by 

F(£) = -Vn + B VQ + Q VB V$V • (p£) 
+ V • [p£v • Vv - pvv ■ V£], 

with the Eulerian perturbation of the total pressure, 



(7) 



(8) 



n = - 7 pV ■ £ - £ • Vp + B • Q, (9) 
and the Eulerian perturbation of the magnetic field, 



Q = V x (£ x B). 



(10) 



By assuming an exponential time-dependence e _lwt for 
the displacement field, we can distinguish two cases. In the 
case of no equilibrium flow, the force operator F(£) is self- 
adjoint, meaning that its eigenvalues u> are purely real or 
imaginary. This results in stable, damped and unstable 
modes. When there is flow, the force operator is no longer 
self-adjoint, meaning that its eigenvalues u> are complex in 
general. This introduces the possiblity of damped stable 
waves as well as overstable modes. 



z.Z. system ot first order ditterential equations 

The Frieman-Rotenberg equation JJJ will be applied in the 
case of a one-dimensional cylindrical plasma equilibrium. 
For this kind of equilibrium, the MHD equations Q— @ 
reduce to the radial force balance equation, 



P9, 



(11) 



where the prime indicates the derivative with respect to 
r. The symbols B, Bg, and vg are the total magnetic field, 
toroidal magnetic field, and toroidal velocity, respectively. 
Furthermore, g represents the gravitational acceleration 
at the distance r, 

GM. 

9 = —> (12) 

with G the gravitational constant and M* the mass of the 
central object. 

For the three-dimensional perturbations, we choose 
Fourier mode solutions of the form 



£z,rnk(r) 



^{mB+kz—ujt) 



(13) 



where m and k are the toroidal and axial wavenumber, 
respectively. This choice can be made because of the sym- 
metry of the equilibrium. Also by exploiting a projection 
based on the magnetic field lines, the Frieman-Rotenberg 
equation Q can be reduced to a system of first order dif- 
ferential equations, which reads 



where 
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In these expressions we use the following definitions: 
k = [Q,m/r,k], k = |k |, the Doppler shifted 



trcqucncy uj = lo — kq ■ v, and b = ko • B. .Notice that 
we have suppressed the subscripts m and k. The 
system (11411 has been derived by Hameiri l|1981J) and 
Bondeson ct al. (1987) for a cylindrical plasma with flow 
but without gravity. Recently, this system was obtained 
by Keppens et al. l(2UU2|> . 

The differential equations l|14|) become singular when 
A = or S = 0, which results in the MHD continua. When 
A = 0, the eigenfrequency u is equal to the local Doppler 
shifted Alfven continuum frequency, 



k 'viwA, 



(21) 



where loa = F '/ J~p is the local Alfven continuum fre- 
quency for static equilibrium. If S = 0, the eigenfrequency 
lo is equal to the local Doppler shifted slow continuum fre- 
quency: 



k -V±LJ S , 



(22) 



where uis = y/ 'jp j '(jp + B 2 ) F/^/p is the local static slow 
continuum frequency. These MHD continua form the ba- 
sic organizing structure of the full MHD spectrum, and 
are schematically shown in Fig. ^ where the equilibrium 
is assumed to be weakly inhomogeneous. The Eulerian en- 
tropy continuum (Goedbloed et al. I2()()4|l fio = ko • v is 
also included in the figure. 



Fig. 1. Continuous parts of the MHD spectrum for a 
weakly inhomogeneous equilibrium. 



2.3. Local dispersion equation 

One can study many of the occuring instabilities by means 
of the local dispersion equation, 



,A 2 S 2 



C + DE 



ASD 



(23) 



where q is the local radial 'wavenumber'. This local dis- 
persion equation can be derived from the differential 
equations (|14|l using the WKB method as discussed in ap- 
pendixE] The equation will be solved neglecting the term 
proportional to ASD jr. This approximation is valid if one 
assumes kc s ^> uj, where c s = \fjpjp is the sound speed. 

The remaining local dispersion equation can be rewrit- 
ten as a sixth-order polynomial in uj which governs all 
discrete local modes. This dispersion equation reads: 



uj 6 + a^uj 4 



■ C13UJ + a-iio + diU) + do = 0, 



(24) 



where the coemcients a, are dchned as 
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with rotation frequency = vg/r, deviation from a 
Keplerian disk V g = Vg/r — g and epicyclic frequency n 2 = 
2ve(rve)' /r 2 . The local dispersion equation reduces 
to the local dispersion equation derived by Dubrulle and 
Knobloch (1993) if one considers axisymmetric perturba- 
tions, a constant density, a zero radial 'wavenumber' q and 
takes the incompressible limit. 



2.4. Magneto-rotational instability 

By making some extra assumptions, we can reduce the lo- 
cal dispersion equation (|23|l further. Consider axisymmet- 
ric perturbations and the situation where qc s , kc s , qv& 
and kvA 2> £ and a purely axial magnetic field. Here, 
v\ = B z j ' J~p is the Alfven speed. In that case, the local 
dispersion equation (|2"4*}> reduces to the following 4 th order 
polynomial: 



& 4 w 4 + b 2 cj 2 + & = 0, 



(30) 



where the coemcients are 
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This 4 th order polynomial can be solved analytically, 
which results in the stability criterion 



El 
p 



> 



rQ 2 ' + r -V g 



P_ 

P 



IP 



(34) 



This is the stability criterion for the axisymmetric MRI 
for disks with a purely axial magnetic field. This crite- 
rion is more generally applicable than the one derived by 
Balbus and Hawley l|1991H to which it reduces for a weakly 
magnetized disk (p> B 2 ). 

3. The accretion disk model in the cylindrical limit 

In order to quantify instabilities using a linear analysis, 
one has to consider the equilibrium state of an accretion 
disk. The model we consider uses power-law scalings for 
the different flow variables, following the self-similar mod- 
els of Spruit et al. I|1987|l . In order to have a model that 
is in an equilibrium state, these profiles have to satisfy 
equation (11). We use the following profiles for the den- 
sity p, thermal pressure p, toroidal magnetic field Bg , axial 
magnetic field B z and the toroidal velocity vg: 

(35) 
(36) 

(37) 

(38) 
(39) 



p = r-V 2 , 
p = c 2 r- 5 / 2 , 
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where the a-parameters express the ratio of toroidal or 
axial to the total magnetic field, 
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and 

V Q 2 = GAh - f + ^ [5(1 +0) (a? + a|) - 4a?] . (42) 

Finally, the parameter /3 measures the radially constant 
ratio between the thermal pressure and the total mag- 
netic pressure , i.e. (3 = 2p/B 2 . Notice that the density 
and pressure profiles are such that the entropy is constant 
throughout the disk, thereby excluding convective modes. 

One of the key parameters of our model is the pa- 
rameter e, which is taken as a constant free parameter. 
For values of e <C 1, the physical interpretation is the 
same as in the model of Shakura & Sunyaev (p[973), i.e. 
e = c s /ve ~ H/r, where H is the scale height of the disk. 
These cases correspond to thin Keplerian rotating accre- 
tion disks. For values of e ~ 1, one is in the regime of 
sub-Keplerian rotating disks (see e.g. Narayan & Yi 1994) 
and for these cases we stick to the interpretation e ~ H/r, 
noting that the identification e = c s /vg is not valid any- 
more. Typically c s /vg will be larger, up to a maximum 
factor of order ~ 10. 

Models of accretion disks which only depend on the 
radius r are sometimes referred to as accretion disks in 
the cylindrical limit (see for example Armitage 119981 and 
Hawley 2001). The model we use is an example of such a 
model, where the results from these models approximate 
the interior equatorial region of the accretion disk. This 
approximation is correct as long as the axial wavenum- 
ber k is much larger than the inverse of the scale height 
H . Therefore the identity k ^> 2n/(er) has to be satisfied 
in order to connect the results from our spectral analysis 
with MHD stability properties of the interior of an accre- 
tion disk. 

4. Numerical codes 

For the investigation of the stability properties two numer- 
ical codes, LEDAFLOW and a Local Dispersion Equation 
Solver (LODES), have been used. In the next two subsec- 
tions we briefly explain their algorithmic details. 

4.1. The spectral code LEDAFLOW 

The code LEDAFLOW developed by Nijboer et al. l(TW7J> 
can compute the complete MHD spectrum for a given 
one-dimensional stationary equilibrium, axial wavenum- 
ber k, and azimuthal wavenumber m. Instead of numeri- 
cally solving the system ()14(l the full set of MHD equations 
Q-lQ} are linearized and then discretized using finite el- 
ements in the inhomogeneous direction, and solved using 
the Galerkin method (Schwarz 1985J)- For the elements, a 
combination of cubic Hermite and quadratic elements for 
the perturbed quantities is used to prevent the creation of 
spurious eigenvalues. The result is a non-Hermitian eigen- 
value problem, which is solved using a QR method or an 
inverse vector iteration method (Press et al. I1988fl . The 
latter one can be used to calculate the eigenfunctions for 
a single eigenvalue. The used boundary conditions treat 



the edges ot the disk as perfect conducting walls. Strictly 
speaking, these boundary conditions are inappropriate for 
considering waves modes in an accretion disk. However, 
interior local modes that do not depend strongly on the 
precise variations in a boundary region are barely affected 
by these conditions. 

4.2. The local dispersion equation solver 

The local dispersion equation solver LODES directly com- 
putes the roots of the local dispersion equation (|2"4*|l at 
a given radial position. This is done by making use of 
Laguerre's method (Press et al. [1988 ). 

The code calculates the local value of coefficients 
a-i H25() - H29|l depending on the given equilibrium, the ra- 
dial wavenumber q, the azimuthal wavenumber m and the 
axial wavenumber k. In principle the code can compute 
all six roots, but there is a switch to calculate only the 
mode with the largest growth rate. Because the code uses 
Laguerre's method (Press et al. I1988|) . it needs an initial 
guess for this root. As an initial guess, we use the analyt- 
ical solution of the most unstable mode of the 4 th order 
polynomial i|30|) . 

4.3. Coupling LEDAFLOW with LODES 

As described in the previous subsection, LODES needs 
a value of the radial 'wavenumber' q and a radial posi- 
tion Ti to calculate the roots. In this subsection, we show 
how we extract these input parameters from the results of 
LEDAFLOW. 

Fig. |U shows the complete MHD spectrum of a calcu- 
lation performed by LEDAFLOW with the following set 
of parameters: m = 0, k — 200, f3 = 1000, e — 0.1 and a 
purely axial magnetic field («i = and eti — 1). In this 
calculation, we use 101 grid points in the radial direction 
on the domain r = [1,2]. It is clear from the resulting 
MHD spectrum that the disk is not stable. We identify 
these unstable modes as MRI modes. Notice that, for a 
purely axial magnetic field, the eigenvalues uj 2 are real 
(i.e. no overstability in this case). 

The eigenfunction \ °f t ne mode with the largest 
growth rate is shown in Fig. From this eigenfunction 
the location r\ and the radial wavenumber q of the in- 
stability can be extracted. The location rj is the radius 
where x has its extremum. The radial wavenumber q fol- 
lows from q = 2n/X r , where A r is chosen as the radial 
wavelength from the inner boundary to the point when the 
value of the eigenfunction is 0.3% of its extremum value. 
The justification of this heuristic method to calculate the 
radial 'wavenumber' q can be found in the appendix, where 
the WKB solutions are matched to the analytical solutions 
in the turning point regions. 

The mode shown in Fig.[3]has also been investigated on 
the enlarged radial domains r = [0.966,2], r = [0.938,2] 
and r = [0.915, 2]. The inner boundary of these domains is 
such that the mode has its extremum at the same radial 
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Fig. 2. Complete MHD spectrum of a weakly magnetized 
thin disk for wavenumbers m = and k = 200. 
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Fig. 3. The eigenfunction x of the mode with the largest 
growth rate. The growth rate is 0.72817 . 



location for the different radial domains. The resulting 
growth rates are all almost equal to each other, which 
confirms our earlier statement that the mode does not 
strongly depend on the boundary conditions used. 



5. The magneto-rotational instability 

We present results on magneto-rotational instabilities in 
an accretion disk in the cylindrical limit using the methods 
described in the previous section. The calculations using 
LEDAFLOW are performed on a radial domain r = [1, 2]. 
We only consider axisymmetric perturbations (m = 0) and 
take a ratio of specific heats 7 = 5/3. We discuss seven 
calculations of which the different sets of parameters are 
shown in Tables and [3 
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5.1. Influence of the toroidal magnetic field 

First, we investigate the influence of the toroidal magnetic 
field on both the growth rate and the oscillation frequency 
of the MRI. This is done by varying the toroidal mag- 
netic field strength, while keeping the axial magnetic field 
strength constant. This means that the plasma beta varies 
between the different calculations. The parameters of the 
four calculations are shown in Table ^ In these calcula- 
tions we consider a thin disk by taking e = 0.1. 

The results are shown in Figs. 14171 In each figure, the 
left panel shows the growth rate and the right panel shows 
the oscillation frequency of the most unstable axisymmet- 
ric mode, both as a function of the axial wavenumber. We 
have scaled the growth rate and the oscillation frequency 
with respect to the rotation frequency. Instead of the ax- 
ial wavenumber k, we use the corresponding local Alfven 
frequency loa scaled with respect to the local rotation fre- 
quency. This scaling is similar to the one used by Narayan 
ct al. (HI!J). 

In our figures, the diamonds represent the calculations 
performed with LEDAFLOW. The dashed and solid line 
correspond with the results from respectively the analyti- 
cal solution of the approximate 4 th order polynomial l|30[) 
and LODES. Recall that the last two calculations make 
use of LEDAFLOW results, in the manner described in 
section ET31 

In the calculations A, B, and C (shown in Figs. 0] 
and El respectively) we change from a purely axial (A) 
to a predominantly toroidal (C) magnetic field configura- 
tion. One sees that the growth rate of the analytical solu- 
tion of 4 th order polynomial (|30|l and LODES matches 
with the LEDAFLOW calculations. The differences to 
LEDAFLOW are less than 4% and 1%, respectively. These 
calculations (A, B and C) also show that the oscillation 
frequency increases away from zero as the toroidal mag- 
netic field strength increases. A similar result was obtained 
by Dubrulle and Knobloch i|1993[l . These authors consid- 
ered an incompressible plasma with a similar toroidal ve- 
locity profile as used in this paper, but a different toroidal 
magnetic field profile. Again, there is a perfect match be- 
tween the LEDAFLOW and LODES results. However, the 
4 th order polynomial (|30|l cannot be used to compute the 
oscillation frequency since it only yields a value for the 
growth rate. The differences between the LEDAFLOW 
and LODES results are again less than 1%. 

Fig. shows results from a calculation for which 
the accretion disk is close to equipartition in a strong 
toroidal \B$\ — 30\B Z \ magnetic field configuration (disk 
model D). Again, there is a perfect match between the re- 
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Fig. 4. The growth rate and the oscillation frequency of the most unstable axisymmetric MM mode as a function of 
the scaled Alfven frequency for disk model A, where (3 — 1000. In this purely axial magnetic field configuration, the 
axisymmetric modes are purely exponentially growing in time. 




Fig. 5. The growth rate and the oscillation frequency of the most unstable axisymmetric MM mode as a function of 
the scaled Alfven frequency for disk model B, where (3 = 500. 



suits from LODES and LEDAFLOW (the discrepancy is 
less than 1%), but the 4 th order polynomial (|30[l no longer 
predicts a correct approximation of the growth rate. This 
is due to the dynamical importance of the toroidal mag- 
netic field component. Comparing Fig. [7] with Fig. El it 
is clear that the oscillation frequency remains roughly the 
same and the growth rate decreases slightly. The analyti- 
cal work performed by Dubrulle and Knobloch (1993 ) also 
indicates a decrease in the growth rate, as a result of the 
increasing tension of the toroidal magnetic field. Hence, 
the cases A-D clearly demonstrate that the MM remains 
active when the toroidal magnetic field increases. This re- 
mains true at least up to equipartition (Fig. [TJt . 

5.2. Keplerian to sub-Keplerian disk 

In this subsection, we investigate the properties of the 
axisymmetric MM in sub-Keplerian disks. Three calcula- 
tions have been performed, all with (3 — 1000, ot\ = 10 
and ct!2 = 1 (see Table EJ. A sub-Keplerian disk can be 
obtained by increasing the value of the e parameter be- 



Table 2. Parameters of Keplerian to sub-Keplerian accre- 
tion disks 



Sim. no. 


e 




Q!2 


P 


k 


Ak 


E 


0.1 


10 


1 


1000 


320-4000 


40 


F 


0.3 


10 


1 


1000 


96-1200 


12 


G 


0.5 


10 


1 


1000 


40-500 


5 



yond e = 0.3. All calculations have been performed such 
that the domain of the Alfven frequency over the rotation 
frequency remains the same: the corresponding domain of 
k- values can be found in Tabled 

The results are shown in Figs. 1811(71 Again, the growth 
rate (left panel) and oscillation frequency (right panel) are 
scaled with respect to the rotation frequency. These scaled 
quantities have been plotted as a function of the ratio of 
the Alfven frequency to the rotation frequency. 

Again, there is a perfect match (less than 1% dis- 
crepancy) for the growth rate obtained with the three 
methods: LEDAFLOW, LODES and the 4 th order polyno- 
mial 1)30(1 . Notice that the scaled growth rates in Figs. 1811111 




LEDAFLOW 
LODES 



Fig. 6. The growth rate and the oscillation frequency of the most unstable axisymmetric MRI mode as a function of 
the scaled Alfven frequency for disk model C, where (3 — 9.901. 
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u A / Q " A / 15 

Fig. 7. The growth rate and the oscillation frequency of the most unstable axisymmetric MRI mode for a disk close 
to equipartition with a dominant toroidal magnetic field component (disk model D). 



coincide although the deviation from a Keplerian disk is 
significant in cases F and G. This can be explained by 
investigating the solution of the 4 th order polynomial in 
the case of a weakly magnetized disk. The solution, scaled 
with respect to the rotation frequency, reads: 

n 2 n 2 2( q 2 + k 2 )\n 2 n 2 y P 9 1P J 

(43) 

For the chosen power-law of the toroidal velocity, the 
epicyclic frequency k 2 = fl 2 . In the case of a weakly mag- 
netized disk and a constant entropy, the term p'V g /p — 
pV 2 /(jp) can be neglected. Furthermore, notice that for 
a pure hydrodynamical case this term is the Brunt- Vaisala 
frequency. For the calculations in this paper, this fre- 
quency is equal to zero because we consider disks with 
constant entropy. 



This results in a equation that only depends on the 
ratio wa/H. Remember that the domain of this ratio has 
been kept the same in calculations E, F and G. 

The oscillation frequencies in Figs. ISIlUl obtained from 
LEDAFLOW and LODES match perfectly (less than 1%). 
These figures show an overall increase in the value of the 
scaled oscillation frequency as one increases the value of 
e. Hence, the overstable nature of the MRI is even more 
significant in sub-Keplerian disks with a dominant toroidal 
magnetic field component. 

6. Conclusions 

We have considered the magnetorotational overstability in 
magnetized accretion disks with a toroidal magnetic field 
component. We have used results from the spectral code 
LEDAFLOW and compared them with more approximate 
solutions of both the fourth and sixth order polynomial lo- 
cal dispersion equation. In our calculations, we have con- 
sidered axisymmetric perturbations (m = 0) within a disk 
in the cylindrical limit. The most important results are 
summarised below: 



ijruwir lULt usliiiu uun requency 

0.81 1 1 1 I 1 1 1 I 1 1 1 I 1 1 1 I 1 1 1 I 1 1 1 I 1 1 1 I I 0.010 1 1 1 I 1 1 1 I 1 1 1 I 1 1 1 I 1 1 r 
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0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 

u A / !) w A / n 



Fig. 8. The growth rate and the oscillation frequency of the most unstable mode for disk model E, where the disk 
rotation is Keplerian. 
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Fig. 9. The growth rate and the oscillation frequency of the most unstable mode for disk model F, where the disk 
rotation is weakly sub-Kcplcrian. 
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Fig. 10. The growth rate and the oscillation frequency of the most unstable mode for disk model G, where the disk 
rotation is strong sub-Keplerian. 



— Magnetorotational instabilities are present in both 
sub-Keplerian and Keplerian rotating accretion disks. 

— The magnetorotational instability is also present in 
Keplerian rotating disks close to equipartition, when 



the toroidal magnetic field strength dominates the to- 
tal magnetic field strength. 
— The oscillation frequency of the instability is nonzero 
due to the presence of a toroidal magnetic field com- 



poncnt tor both sub-Kepierian and Keplerian rotating 
disks. 

In an accompanying paper by van der Swaluw et al. 
we will consider accretion disk models that also 
allow for convective instabilities, and analyse their relation 
with respect to MRI modes. All our solvers can be used 
to study non-axisymmetric perturbations (m 7^ 0). This 
is left to future work. 
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Fig. A.l. The radial 'wavenumber' q of the eigenfunc- 
tion x shown in Fig. EI The growth rate is 0.72851 . 



of ordering yields the following expressions for the radial 
'wavenumber' q and the amplitude p: 




1 



(/s) 1/4 ' 



(A.6) 



Appendix A: Radial 'wavenumber' 

To justify how we calculate the radial 'wavenumber' q as 
discussed in subsection we discuss the WKB approxi- 
mation in detail. Instead of using the system of first order 
differential equations (fH|) we use the equivalent general- 
ized Hain-Lust equation, 



{fx')' + 9X = 0, 



where 



= AS 
1 ~ rD" 



ASD 



C + DE 



ASD (C_ 
D 



(A.l) 

(A.2) 
(A.3) 



This equivalent equation is preferred over system (|14fl for 
a more straightforward WKB analysis. In standard WKB 
analysis for the generalized Hain-Lust equation (|A. 1|) we 
write the solution in the following form: 



Using the expression for q, we can plot this quantity 
for a given radial domain. This is shown in Fig. IA.1I for 
the eigenfunction \ plotted in Fig. [31 Here, we applied the 
same radial domain as used in the LEDAFLOW calcula- 
tion. The plot shows that at a certain radius r\ the radial 
'wavenumber' becomes zero. Close to this radius the WKB 
approximation fails, because qL ^> 1 is no longer satisfied. 
However, we can find an analytical solution close to this 
radius. To find the analytical solution we write the gen- 
eralized Hain-Lust equation (|A. 1|) in the following form 



'X 



dR 2 



G X = 0, 



where 



R'=J, and G = fg. 



(A.7) 



(A.; 



Both the radial 'wavenumber' q and the function G be- 
come zero at the radius r\ as expected. We apply a Taylor 
expansion of G about Ri(= R(ri)) and insert this in the 
generalized Hain-Lust equation (|A.7J| . This yields 



X(r) = p{r) exp 



q(r)dr 



(A.4) 



fx 

dz 2 



-— - ZX = 0, 



(A.9) 



where tq is the radial location of the inner boundary. 
Inserting this solution in the generalized Hain-Lust equa- 
tion (lA.lll gives 

-fq 2 p + 9P + fp" + f'p' + i (2JV + f'qp + fq'p) = 0. 



where 



dG 



dR 



1/3 



(R - i?r) . 



R=R t 



(A.10) 



(A.5) 

Airy function, 

Introducing a length scale L for the radial background 
variation, we can neglect the second order terms fp" + f'p' 
compared to the first order ones if qL ^S> 1. This kind 



The analytical solution of this differential equation is the 



Ai(z) = - 



cos (|s 3 



sz) ds, 



(ATI) 



which tor large \z\ has the asymptotic torm 
1 



Ai(z) 
Ai{z) 



2^z 1 / i 
1 



exp 



.2 3/2 



We now calculate the integral 



" q {r)dr= 2 -{-zf'^ 



which, for R < 0, results in the solution 



X(r) 



An 



{fa) 1 



q(r)dr + — tt 



for z > 0, 

for z < 0. 
(A.12) 

(A.13) 
(A.14) 



where ^4o is constant. We apply to this function the 
same boundary condition as in LEDAFLOW at the in- 
ner boundary. This means that x( r ) is zero at the inner 
boundary, which yields 



q(r)dr 



1 



■ 717T, 



(A.15) 



where n is a natural number. If we numerically calcu- 
late this integral for the radial 'wavenumber' q plotted 
in Fig. IA.1I we find that it is approximately 0.767T. For 
numerical integration we used Simpson's rule (Press et 
al. I1988JI . Thus determining the radial 'wavenumber' as 
discussed in subsection 14.31 gives approximately the right 
value for q. Hence, even for this most unstable mode, the 
local WKB analysis is in excellent agreement with the nu- 
merical solution. 
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